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The inverse Ising problem consists in inferring the coupling constants of an Ising model 
given the correlation matrix. The fastest methods for solving this problem are based on 
mean-field approximations, but which one performs better in the general case is still not 
completely clear. In the first part of this work, I summarize the formulas for several mean- 
field approximations and I derive new analytical expressions for the Bethe approximation, 
which allow to solve the inverse Ising problem without running the Susceptibility Propagation 
algorithm (thus avoiding the lack of convergence) . In the second part, I compare the accuracy 
of different mean field approximations on several models (diluted ferromagnets and spin 
glasses) defined on random graphs and regular lattices, showing which one is in general more 
effective. A simple improvement over these approximations is proposed. Also a fundamental 
limitation is found in using methods based on TAP and Bethe approximations in presence 
of an external field. 

Mean-field approximations (MFA) are very important tools in statistical mechanics, since they 
provide an approximated description of a physical system in terms of few parameters (e.g. local 
magnetizations). Among MFA the one based on the Bethe approximation (BA) is very effective. 
In recent years the BA — originally derived for the ferromagnetic model on regular lattices [l| - 
has been extended, under the name of Cavity Method, to models having arbitrary couplings and 
topologies Although the BA is exact only for tree-like topologies, its application to models 
defined on random graphs has proved very successful: see e.g. the cases of low-density parity check 
codes, spin glasses and constraint satisfaction problems, all nicely reviewed in Ref. y. 

The inverse Ising problem, originally known as Boltzmann machine learning, consists in inferring 
coupling constants of an Ising model (both pairwise interactions and external fields) given the vector 
of magnetizations and the matrix of pairwise correlations. In recent years the inverse Ising problem 
has received a lot of attention, specially in connection to inference in biological problems 4tlL\- 

The inverse Ising problems can be viewed as the dual problem with respect to the 'direct' 
problem of estimating magnetizations and correlation given the Hamiltonian. So, under any MFA, 
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the inverse problem can be solved by inverting (if possible) the analytic expressions that give 
the magnetizations and the correlations as a function of interactions and fields. Although MFA 
usually do not provide directly the correlations between distant variables, these correlations can 
be computed by using the linear response theorem ^]. 

Within the BA, a very fast and efficient way to estimate correlations between any pair of 
variables is given by the Susceptibility Propagation (SuscProp) algorithm recently introduced in 
Ref. 0. SuscProp is an iterative algorithm for solving a set of self-consistency equations: when it 
converges is very fast, but sometimes it may not converge. Indeed the application of the BA to the 
inverse Ising problem has been limited up to now by the range of convergence of SuscProp 
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In this work I present an analytical expression for the fixed point of SuscProp, thus avoiding 
any problem related to its lack of convergence. Actually such an expression already appeared in 



Ref. 



12j, but was unknown to the statistical mechanics community (including the author): otherwise 



there would be no need for the SuscProp algorithm introduced in Ref. |9|. 

In the first part of this work I derive new analytical expressions for solving the inverse Ising 
problem under the BA. These analytical expressions, allow for a fair comparison among different 
MFA, in a wide range of temperatures, both for the problem of estimating 2-point correlations 
given the couplings (direct problem) and for the problem of estimating couplings given correlations 
and magnetizations (inverse problem). 
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20j ] and not only the correlation 



In the present work I consider MFA obtained from the so-called Plefka expansion 
from a small correlations expansion [16]. For the inverse Ising problem I compare methods that 
take in input only the correlation matrix. More complex, but usually slower, inference methods 
exist that require many samples of equilibrium configurations 
matrix. 

Improving over these MFA would be very welcome. It is well known that MFA ignore loops, 
so correcting these MFA by adding the loops contributions would be the right direction to follow, 
rlowever at present all the methods which has been developed to consider explicitly the loops 
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251 ] do not provide analytical expressions for the correlations, which are simple enough to be 
inverted. For this reason, I have not considered these improved algorithms in the comparison of 
MFA for solving the inverse Ising problem. 

Nonetheless, I am proposing a simple improvement of inference methods (both for the direct 
and the inverse problems) based on the idea that the loops may modify similarly self-correlations 
and correlations between close-by variables. 
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THE MODEL AND THE MEAN-FIELD APPROXIMATIONS 



In order to keep the presentation simple, I prefer to deal only with binary variables (Ising 
spins) Si = ±1 and Hamiltonian containing up to two-body interactions, i.e. external fields and 
pairwise couplings. Thus, the most general model I want to study is defined by the following joint 
probability distribution over N Ising variables 



P(si, 



1 



,SN) 



Z(J,h) 



exp 
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(1) 



where the partition function Z(J,h) is a normalizing constant, that depends on all the couplings 
J = {Ji j} and the external fields h = {hi}. Please notice that the temperature parameter has 
been absorbed in the definition of external fields and couplings. All the required information about 
the model is encoded in the free-energy 



F(J,h) = hxZ(J,h) . 



(2) 



In the rest of this Section I summarize the most common MFA to the free-energy: I am particularly 
interested in deriving the self-consistency equations for the magnetizations that are used in Section 
|TI]for obtaining 2-point correlations. 

The simplest MFA, also known as naive MF (nMF), approximates the model in terms of local 
magnetizations rrii = (sj), where the angular brackets represent the average w.r.t. the measure in 
Eq. ([!]). The corresponding approximation to the free-energy is 



-FnMF = 
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where H(x) = — xln(x) and the mi must be fixed according to the self-consistency equations 



(3) 



dm; 



Jijrrij + hi — atanh(mj) =0 => mj = tanh 



A better MFA can be obtained by considering also the Onsager reaction term 
the following TAP approximated free-energy and self-consistency equations 
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rrii = tanh 



hj + Jijirrij - Jij{l - m|)m, 



(5) 
(6) 
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In the TAP approximation, when computing the marginal probability of spin Sj (i.e. its mag- 
netization rrii), the reaction term modifies the marginal probabilities of the neighboring spins, 



mi 



(rrij - Jjj(l - m|; 



has been recognized 
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), in order to try to remove the effect of the spin Sj under study. It 



14] that F„ MF and F XAP are only the flrS t two te„„ S of the expansion 
of F(J,h) in small couplings J at fixed magnetizations m = {rrii}. This expansion contains [14| 
both loop terms, like JijJjtJu, and terms with higher powers of a single coupling, i.e. J^: the 
latter terms, that correspond to considering recursively the reaction to the reaction between spins 
Si and Sj, can be resummed and lead to the BA. 

The BA gives a description of the model in terms of magnetizations rrii and connected correla- 
tions Cij = (siSj) — rriimj between neighboring spins (i.e. spins connected by a non-zero coupling 
Jij). The BA can be derived in two equivalent ways. The first way consists in finding values of m 
and c minimizing the following free-energy 

(1 - m*)(l - rrij) + c^' 



H | {I + mi){l + rrij) + Cij +R 



+ H 



(1 + rrii)(l - rrij) - dj 
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where cij is the degree of spin Sj, i.e. the number of its neighboring spins. In Eq.Q the last two 
terms correspond to the average value of the energy at given magnetizations and neighbouring 
correlations, while the first two terms correspond to the entropy of the Bethe approximation to the 
joint probability distribution of the spin variables, 

BA -n- pjj(sj,Sj) 



P(si, 



n 

(ij) 



Pi{Si)pj{Sj 



(8) 



where the first product runs over all pair of neighboring spins and the two-spins and single-spin 
marginal probabilities are given respectively by Pij(si, Sj) = [(1 + mjSj)(l + rrijSj) + CijSiSj]/4 and 
Pi(si) = (1 + rriiSi)/2. The conditions dF-Q^/dcij = can be solved analytically and lead to 



Cij {rrii , rrij , tij ) 




((1 + mj)(l + rrij) + Cij \ {{l - mj)(l - rrij) + CiA 



(1 + mj)(l - rrij) 

2 



(1 - rrii)(l + rrij) - Ci 



(9) 
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y(l~ ^ij) 2 ~ ^Ujim - tijmj)(rrij - - wijmj . (10) 



where = tanh(Jy). Please note that Eq.Q is identical to Eq.(26) in Ref. 



16 and this is a 



further confirmation that resumming all 2-spin terms in the Plefka expansion leads to the BA. 



Moreover Eq.([9]) has been used in the literature [3, [3] as the independent-pair (IP) approximation 
for inferring couplings from magnetizations and correlations: such an approximation infers the 
coupling Jij by assuming spins and Sj form an isolated pair with magnetizations m; and rrij and 
correlation c«. Unfortunately under this IP approximation computing the external fields in not 
immediate and moreover even the estimates of the couplings are rather poor (see Section |V|) , 

By making the substitution Cy — > (m, , rrij , tij ) in Fba one can obtain the Bethe free-energy 
only in terms of magnetizations, from which the self-consistency equations for the magnetizations 
can be derived. However this derivation requires a rather complicated algebra and I prefer to 
obtain the same equations in a much simpler alternative way. 

In the so-called Cavity Method [2] local magnetizations rrij and neighbouring correlations Cjj 
are expressed in terms of some auxiliary variables, the cavity magnetizations m^ (i.e. the mean 
value of Si in the absence of a neighboring spin Sj): 

m = 77 W+ 77 ' 

1 + m>(' m>- 

1 ' ' (12) 



1 , {3)4. (0 ' 
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* , U) (*) 
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Cavity magnetizations must satisfy the self-consistency equations 



(i) 

m) = tanh 



hi + atanh(t ifc rn£ 



(0> 



(14) 



28]: 



These equations are often solved by an iterative algorithm known as Belief Propagation (BP) 
in case of convergence, the fixed point of BP gives directly the Bethe free-energy that admits an 
expression in terms of cavity magnetizations only Q]. 

In order to obtain a closed set of self-consistency equations in the magnetizations m, let me 
solve eqs. (jimi2|) for the cavity magnetizations and find 

my = f (mi, rrij, Uj) my = f(mj,mi,tij) , (15) 

where 



1 - t 2 - y/(l - t 2 ) 2 - 4t(mi - m 2 t)(m 2 - mTt) 
/(mi ' m2 ' t] = 2t(m 2 - mTJ) • (16) 

The sign in front of the square root has been chosen such that f(0, 0, t) = as it should. A 
consistency check can be made by substituting expressions (fT5j) in Eq. lfTBj) to obtain again the result 
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in Eg. (jlOp . Finally, combining Eq. ljlip and Eq. (|14p . it is possible to obtain the self consistency 
equation for the magnetizations under the BA: 



rrij = tanh 



(17) 



j 

It is fair to comment that the use of this formula for finding Bethe magnetizations is not a good 
idea: indeed an iterative solution of Eq. (|17|) is typically more unstable than BP solving Eq, (|14|) . 
My interest in this formula is that it involves only physical magnetizations (not cavity ones) and 
can be used to obtain correlations (see Section [XT]) and to solve in a fast way the inverse Ising 
problem (see Section |V|) , 

A series expansion of the exponent in Eq. (|17p for small couplings gives 

h + ^ a,tanh(tijf(mj,mi, tjj)j ~ hj + ^ (jjjmj - J?-(l - m^rm + ■ ■ ■ ) , (18) 
j j 

and one recognizes that the first two terms of the expansion are the naive MF approximation and 

the Onsager reaction term. This expansion should make clearer that the BA is a way of considering 

recursively all the reactions between a pair of neighboring variables. 



II. COMPUTING CORRELATIONS BY LINEAR RESPONSE 

A preliminary step to solve the inverse Ising problem by any MFA is to derive an analytical 
expression for the pairwise correlations as a function of the coupling constants. Actually, the MFA 
discussed in Section [I] do not provide information about the correlation between distant variables: 
indeed, naive MF and TAP approximations give cj. = for any pair of variables, and the BA only 
provides an expression for correlation between neighboring spins, see Eq. (|10p . which is trivially 
Cij = tij in case of null magnetizations. 

Nonetheless, a closed set of equations for the connected correlations 1 , Cij = {siSj} — (si)(sj) for 
any pair i,j, can be derived from the magnetizations self-consistency equations, Eqs.(HJ), (0, (|17j) . 
through the linear response [3] 



1 Please do not confuse the correlation dj with the parameter cy appearing in the BA: the two coincide only when 
the BA is exact. 
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The inverse correlation matrices C 1 for the three MFA discussed above are given by the following 
expressions: 



naive MF {C~^ ¥ \ 
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(20) 
(21) 

(22) 



where fi(mx,m2,t) = df(mi,m2,t)/dmi and f2( m i, rn >2,t) = df(mi,m2,t)/dm2- From these 
expressions one can obtain directly any correlation by simply computing the inverse of a matrix. 

Please note that Eq, (|22p gives exactly the same solution found by the SuscProp iterative al- 
gorithm [9], which is presently considered one among the best inference algorithms. The main 
advantage of Eq. (|22p is that it always provides the correlation matrix, even in those cases where 
SuscProp does not converge to the fixed point. Moreover inverting a matrix takes roughly the same 
time of a single iteration of SuscProp, and so using Eq. (|22p is much faster than running SuscProp, 
even when the latter converges. 

Nevertheless, it is fair to notice that the use of Eq. (|22p does not solve all the problems related to 
the lack of convergence of SuscProp. Indeed, during the many tests I have run, I noticed that often 
the lack of convergence of SuscProp does correspond to the BA fixed point becoming unphysical: 
in these cases, by inverting the correlation matrix provided by Eq. (|22p . one gets an unphysical 
correlation matrix (e.g. a correlation matrix with negative diagonal elements!). In this sense the 
lack of convergence of SuscProp gives a warning that the "blind" use of Eq. (|22|) does not provide. 
So, a general suggestion when using the above formulas, providing an analytical expression for the 
correlation matrices under a MFA, is to check explicitly the physical consistency of the outcome. 

One may comment that Eq. (|22|) contains the magnetizations and the iterative computation of 
these (i.e. the BP algorithm) suffers the same convergence problems of SuscProp: this is easy to 
prove, given that the homogeneous SuscProp equations are nothing but the iterative equations for 
evolving under BP a small perturbation in the magnetization, and so BP is unstable if SuscProp 
does not converge. However there are provably convergent algorithms for the computation of 
magnetizations under the BA 
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30 ]: the use of these algorithms in conjunction with Eq. (|22p 
allows a direct computation of correlations under the BA. Moreover there are situations where 
magnetizations are known a priori and Eq. (|22p can be applied directly: e.g. when symmetries in 
the probability measure force magnetizations to be zero, or in the inverse Ising problem, where 
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magnetizations are given as an input to the problem. In the rest of the paper I deal mainly with 
these two cases. 



A. Estimating correlations in case of null magnetizations 



A preliminary ranking of MFA can be done on the basis of how good are their estimates of 
correlations, given the couplings. Indeed I expect that the better is this estimate, the better will 
be the solution to the inverse problem. 

For simplicity I concentrate on models with no external fields and the couplings are multiplied by 
a parameter f3 (the inverse temperature) such that the difficulty of the inference problem increases 
with /?. 

In case of null magnetizations, the expressions for the inverse correlation matrices simplify a lot 
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since /i(0,0,t) = 1/(1 - t 2 ) and / 2 (0,0,t) = -t/(l -t 2 ). 

Given that for m{ = the expressions for the correlation matrices are much simpler, I report 
also those that can be obtained from the Plefka expansion at the third and fourth order 
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4 th order 
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(27) 



The purpose is to understand if and how much does the estimate of the correlation matrix improve 
by adding terms in the Plefka expansion. 

I have tested the accuracy of formulas in Eqs. ()23ll27|) for ferromagnetic (Jy = 1) and spin glass 
(Jjj = ±1) models defined on fully connected (FC) topologies, on a 2D square lattice and on a 
3D cubic lattice. In diluted versions of these models a fraction (1 — p) of couplings has been set 
to zero. In models defined on FC graphs the couplings have been normalized such as to have a 
critical inverse temperature f3 c = 1 in the thermodynamic limit. 
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FIG. 1: Error made by 5 mean-field approximations in estimating the correlation matrix, given the couplings. 
Shown are typical samples of size N = 5 2 (the qualitative behavior does not change for larger sizes). 



The discrepancy between true correlations C and those inferred C is defined as 



(28) 



i ,3 



In Figure [T] and [2] I report the typical behavior of the error Ac between exact and estimated 
correlation matrices for 5 different MFA. Figure [T] shows results for models defined on a 2D square 
lattice, while Figure [2] refers to FC and 3D topologies. In order to compare the MFA estimates 
with the exact correlation matrices I am studying here small systems, but the qualitative behavior 
does not change for larger sizes. 

Although the quantitative behavior of Ac depends on the specific sample, some general state- 
ments can be made: 

• naive MF is typically the worst MFA and shows many spurious singularities (roughly one 
for each peak in Ac); 

• TAP and 4 th order approximations typically show no (or very rare) singularities; 

• the best estimate is typically provided by BA and TAP, with BA being the best unless it 
has a singularity (in this case TAP becomes the best at lower temperatures, higher /3). 
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FIG. 2: Same as in FigureQ]for typical samples of the fully-connected spin glass (SK model) of size N = 20 
and of the 3D diluted ferromagnet of size N = 3 5 . 

These results suggest that increasing the number of terms in the Plefka expansion does not always 
improve the estimate of the correlation matrix (as one could have naively expected). On the basis 
of these preliminary results I believe it is relevant to consider only TAP and BA for the inverse 
Ising problem, together with other inference methods (see Section ITV)) . 

In the left panel of Figure the results obtained by TAP and BA are almost perfectly super- 
imposed (indeed the former in not well visible). This is expected since the TAP approximation is 
exact for the SK model at high temperatures (/3 < f3 c = 1) in the large N limit: so the BA can 
not improve it, but in 1/N corrections. Indeed a careful analysis shows a tiny improvement of BA 
over TAP around the critical temperature, where 1/N corrections are stronger. 

Please note that in the right panel of Figure [2] the high temperature (small (3) behavior of 
Ac is very different than in previous plots: indeed for j3 — > 0, Ac goes to a constant, instead of 
decreasing with a power law in f3 (as in Figure Q] and in the left panel of Figure [2]). This is due 
to the fact that the comparison has not been made with the exact correlation matrix, but with 
correlations measured from a Monte Carlo (MC) simulation. Actually, in this case, I have used the 
Wolff algorithm and the correlation matrix has been computed from 10 5 independent measures. 
The difference between the error due to the MFA and the error due to MC noisy data can be 
better appreciated in Figure in the high temperature region the error does not decrease below 
a limiting value given roughly by the inverse of the square root of the number of measures. 
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III. IMPROVING INFERENCE ALGORITHMS 



Expressions in Eqs. (|23ti27|) are intrinsically approximated, and turn out to be correct only 
in some particular cases. Naive MF and TAP approximations (as well as 3 rd and 4 th orders 
approximations), being the first orders in a small couplings expansion, are exact only in the limit 
of very weak couplings (either high temperature or fully-connected models in the large iV limit). 
The BA, on the contrary, is exact also for coupling intensities O(l), but only if the interacting 
network is a tree; on random graph models (which are locally tree-like) the BA turns out to be 
correct as long as the model has only one state (modulo the known symmetries). On any other 
model those expressions are approximated and it is worth trying to improve it. 

Let me first notice that any of the above MFA returns in general a value for the self-correlation 
differing from the exact one, i.e. Cu ^ 1 (for simplicity I consider the case of null magnetizations, 
but the argument is general). This fact can be easily explained, noticing that all the above MFA 
assume that correlations along loops are vanishingly small (at least in the large N limit). On the 
contrary, on any loopy graph, like e.g. regular lattices, correlations along loops are important and 
may alter significantly the mean- field estimates. A general solution to this problem is still not 
available, although many work is in progress to include loop contributions to MFA 
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25]. 



What I am proposing here is a simple heuristic improvement. Once the correlation matrix C^ is 
computed by one of the approximations described in Section HH a properly normalized correlation 
matrix can be defined 

Cii - -== ■ ( 29 ) 

y ^ii^jj 

By definition Cu = 1 and also off-diagonal element may approximate better the true correlations. 
The reason for this is that the loops neglected in MFA actually modify in a similar way both 
self-correlations Cu and off-diagonal correlations Cij, and the heuristic normalization in Eq. (j29p is 
assuming that the modifying factor only depends on the loop structure around sites i and j (which 
is certainly wrong for distant sites, but may be a reasonable approximation for closed- by sites). 

In Figure [3] full points show that the error Ac in the BA decreases by roughly one order of 
magnitude if normalized correlations are used. On the contrary, the TAP result is not very sensitive 
to this normalization: the reason is that the estimates of the self-correlations in TAP remain quite 
close to the right value, specially if compared to BA estimates that diverge at the singularity 
(mark by a peak in Figure [3]) . On the right of such a peak the error obtained by the normalized 
BA is not reported because Eq.Q29l) can not be used, since several BA estimates for self-correlations 
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FIG. 3: Same as in Figure Q] for a typical sample of the 2D diluted ferromagnet of size N = 5 2 . The error 
Ac has been computed with respect to the exact correlation matrix and with respect to the one measured 
in MC simulations. Full points show the error obtained with the normalization trick. 

are negative. This is the problem of the BA fixed point becoming (strongly) unphysical, already 
discussed in Section [XTJ indeed by running SuscProp on this sample one would observe convergence 
only for (3 smaller than the peak location. I would like to stress again that checking the physical 
consistency of a solution based on a MFA is very important: for the sample shown in Figure El 
even without knowing the exact correlations, one should switch from the BA to TAP, when the 
former reaches the singularity (that manifests e.g. in SuscProp not converging or in self-correlations 
diverging) 2 . 

Moreover there are cases (e.g. homogeneous FC models) where the spurious singularity induced 
by the MFA in a system of finite size is such that Ca and CV,- diverge with the same law at the 
spurious critical point, while the normalized correlation Cy stays finite (and much closer to the true 
one). For example for the FC ferromagnetic model the normalized correlation Cmfa estimates the 
true correlation with an error roughly half than the one of Cmfa for any of the 5 MFA considered 
here. 

2 Actually for a ferromagnet one knows how to break the up-down symmetry and let BP converge even at low 
temperatures: once BP returns non-zero magnetizations mi, the correlation matrix can be computed by mean of 
Eq. (|22|l . However in the general case, BP does not converge in presence of long range correlations, i.e. after the 
singularity, and one must resort to other MFA. 
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IV. METHODS FOR THE INVERSE ISING PROBLEM 

I consider 4 different approximations for solving the inverse Ising problem. The simplest one 
is the independent-pair (IP) approximation, already discussed in Section U and recalled here for 
convenience 



4f 7 '"I ¥ ~ ,.-{); : 4 I • (30) 
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((1 + mi)(l + rrij) + ((1 - mj)(l - m,) + C^) 



4 



(1 + mj)(l - rrij) - dj\ ({1 - mj)(l + mj) - Cij 



Among the MFA which can be derived from the Plefka expansion, I consider only TAP and BA, 
because are those performing better in the direct problem of estimating correlations (see Section HI]) . 
The corresponding expressions for the inferred couplings can be obtained by solving the equation 

2m i m j Jf j + J ij + (C- 1 )i j = V(t^i) (31) 

for TAP and the equation 

(C-%= ~% f \ ( ^' mu % = . " % V(i^i) (32) 

i - t {j j (m^m^tij) ^ _ j2.) 2 _ 4 tij (mi - /„////)!///, - t ijmi ) 

for the BA, thus leading to 
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2(C- 1 ) i 
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The fourth approximation I am considering has been obtained from a small correlation expansion 
by Sessak and Monasson [Iff] and has been further simplified in Ref. 



271 to the following expression 



J|M_ (C-% + J* ^— . (35) 



For each approximation, I measure the error in inferred couplings J-j with respect to the true 
ones Jij by the following expression 



/-~/i<l in 



(36) 



I study both diluted ferromagnetic model with a fraction p of non-zero couplings (Jij = f3) and 
undiluted spin glass models (Jij = ±/3 with probability 1/2). I also consider several topologies: 2D 
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square lattices, 3D cubic lattices, random regular graphs with fixed degree c = 4 and fully connected 
(FC) graphs. In the latter case the couplings are normalized in order to have a phase transition at 
f3 c = 1 in the thermodynamic limit. I restrict the study to models of small sizes, with N ranging 
between 20 and 100, because these are the sizes for problems of biological interest. Moreover, as 
discusses below, the number M of independent measurements of the correlation matrix that make 
inferred coupling reasonably good grows linearly with N and so for larger systems the number of 
measurements needed become too large. The data shown in Section [V] have been obtained with 
M = 10 6 independent measures of the correlation matrix (unless differently stated) and going to 
much larger values seems to me rather unrealistic if compared with practical applications. 

A. Normalization trick for the inverse Ising problem 

The trick of normalizing the correlation matrix to improve inference (see Section [TTTJ) can be 
extended to the inverse Ising problem. In practice, it corresponds to solve all the equations relating 
the inverse correlation matrix C _1 to the couplings Jij, including also those for the diagonal 
elements which are usually ignored. 

Let me illustrate the new method for the simple case of the TAP approximation with null 
magnetizations. In this case, solving the inverse Ising problem only on the off-diagonal elements is 
equivalent to solve the equations 

but the diagonal equations are in general unsatisfied 

k 

where D is the inverse correlation matrix estimated by TAP once coupling Jij are given. Please 
notice that diagonal elements Da are fully determined once the off-diagonal elements are known. 

Following Eq. (j29p I would like to normalize the estimated matrix D, and produce a normalized 
inverse correlation matrix D, matching better the true inverse correlation matrix C" 1 . In other 
words I would like to solve the equations 

(C-% = = DijXiXj , (37) 

where the N variables Aj are exactly those necessary to solve the new N diagonal equations. 
Physically speaking, A^ should be the self-correlation Ca produced by the MFA, when using the 
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right couplings Jy. When the Aj becomes very different from 1, then the MFA is working very badly; 
however I expect situations where some of the errors produced by the MFA can be compensated 
by this normalization trick. 

Unfortunately the solution to the new equations, those involving both Jj,- and \, does not have 
an analytical expression and need to be solved numerically. I have adopted an iterative solution, 
which is very fast (when it converges). In practice, I start with all Aj = 1 and then, iteratively, first 
I solve the off-diagonal equations, thus getting an estimate for the couplings, and then I solve the 
diagonal equations to obtain a new estimate for the A's. I repeat this iterative procedure, updating 
A's with a damping factor, until the A's variations is below a threshold (typically 10 -8 ). 

From the many tests I have run, I noticed that this normalization trick is more relevant for 
unfrustrated models (as the diluted ferromagnets studied below) or models containing regions very 
weakly frustrated (those usually leading to Griffith singularities). Most probably in these weakly 
frustrated regions correlations get self-reinforced by the loops (ignored in the MFA) and thus 
the normalization trick may improve coupling estimates. On the contrary for strongly disordered 
models, like spin glasses, the effect of the normalization trick depends a lot on the disorderd sample 
and it does not seem to give a clear improvement on average. 

The use of the normalization trick for improving inference in the inverse Ising problem may 
resemble the diagonal-weight trick introduced in Ref. [8|, but it is actually very different. In the 
diagonal-weight trick, the self-couplings Ja are allow to take non-zero values in order to solve all 



the equations (C 1 
It has been shown 



D{j, while in the normalization trick the self-couplings Ja remain null. 



151 ] that the the first order approximation (nMF) with the diagonal-weight 
trick improves over the second order approximation (TAP) in estimating magnetizations. However 
the estimates for the couplings do not improve at all, because the off-diagonal equations are left 
unchanged by the diagonal-weight trick. On the contrary, the normalization trick used here does 
change the estimates for the couplings. Moreover, being based on the very general requirement 
that self-correlations must take the right value, it can be applied to any approximation. 

V. NUMERICAL RESULTS ON THE INVERSE ISING PROBLEM 

Let me start by making some general statements that summarize the numerical results shown 
in this Section. Among the four approximations studied (IP, TAP, SM and BA) it seems in general 
true that: 

• IP always provides the worst estimate, especially at high temperatures (low /3); 



16 



• BA always outperforms TAP; 

• between SM and BA, the former is better in the high temperature (low f3) phase, while the 
latter is better at lower temperatures (higher /3 values); 

• at low /3, the error Aj in couplings inference is completely dominated by the uncertainty on 
the correlation matrix and it does not depend on the inference method: for this reason the 
range where SM is the best method becomes tiny, especially for noisy data; 

• at high /3, the error Aj produced by TAP and SM diverge, the one by IP stays limited, but 
very high, and only BA may have a reasonable error; 

• for diluted ferromagnetic models the normalization trick works fine and thus BA with the 
normalization trick is the best method overall; 

• in presence of an external field, i.e. when magnetizations are different from zero, TAP and 
BA stop working at high enough (3 values (i.e. these methods do not admit a solution); given 
that at high f3 neither IP nor SM provide acceptable inferred couplings, I conclude that in 
such a situation the inverse Ising problem needs to be solved by other methods not explored 
in the present work. 



Let me start discussing the high temperature (low f3) regime. In Figure H] are reported the 
errors Aj in inferring the couplings of a ferromagnetic model of N = 20 variables on a regular 
random graph of fixed degree c = 4. For simplicity I am plotting only the results obtained with 
SM and BA. Data in the upper curves have been computed using a correlation matrix averaged 
over M independent measurements, while data in lower curves have been computed from the exact 
correlation matrix. It is clear that the low j3 regime is completely dominated by the uncertainty in 
the correlation matrix (as already noticed in Ref. [10(]), and the error in this regime is independent 
on the inference method used. For this reason the good performances of SM in this regime are 
actually washed out and, even for M = 10 6 measurements, the improvement of SM over BA is very 
limited (see Figure HJ). Moreover such an improvement tends to become smaller by increasing the 
system size because in the low (3 regime the error goes like 



A. Ferromagnetic models 




(38) 
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FIG. 4: Dependence of the error Aj in inferred couplings on the number M of independent measures for 
the correlation matrix. The model is a ferromagnet of N = 20 variables on a regular random graph of fixed 
degree c — 4, whose critical temperature is marked by the vertical dotted line. 



10 



0.01 



0.1 r, 



ferromagnet N=20 
Bethe lattice (c=4) 


* X n 


v a 

\ * S *' ' 

>fL/ . 


/ 

IP — i ' 

TAP 

SM ••••« 

BA □ 
TAP norm 
BA norm 



10 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

P 



0.01 









ferromagnet N=100 
Bethe lattice (c=4) 


X 

f 


* xj< ,"° D 

□ 

□ 

□ 




; ; s 
: X ■ 

// 








\ / x 

Y 

\ / 4 
■, 


i ; 

, A 
a - 


IP —I — • 

TAP 
SM 

BA □ 
TAP norm 
BA norm - 


3 0.1 0.2 0.3 


0.4 0.5 

P 


0.6 0.7 0.8 0.9 



FIG. 5: Errors in the couplings inferred by several approximations. The model is a ferromagnet on a random 
regular graph of fixed degree 4. Shown are two typical samples of sizes N — 20 (left) and N — 100 (right). 
The vertical dotted lines mark the locus of the ferromagnetic phase transition in the thermodynamic limit. 



I think that comparing inference methods by using the exact correlation matrix is rather un- 
realistic, given that in any practical application the correlations are always known with some 
uncertainty. So in presenting below the numerical results I always consider the case with M = 10 6 
independent measurements for the magnetizations and the correlations. 

In Figure [5] I am showing the error in the couplings inferred by several approximation for a 
ferromagnet on a random regular graph with fixed degree c = 4. The two panels correspond to 
sizes N = 20 (left) and N = 100 (right) and show that the qualitative behavior is mostly size 
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FIG. 6: Errors in the couplings inferred by several approximations. The model is a diluted ferromagnet on 
a 2D square lattice of size N — 7 2 . Shown are typical samples for 4 different dilutions: 0.6, 0.7, 0.8 and 0.9. 
The vertical dotted lines mark the loci of the ferromagnetic phase transitions in the thermodynamic limit. 



independent. Also the dependence on the specific sample (i.e., on the random graph) is rather 
weak. The data in Figure [5] support many of the statements written above: (i) IP is a very bad 
approximation even in the low (3 regime; (ii) BA always outperforms TAP; (iii) SM is better than 
BA only in the low j3 regime, but here the error is dominated by the uncertainty in the correlation 
and increases with the system size, so the improvement of SM over BA is tiny; (iv) errors in TAP 
and SM diverge for large f3, while those in IP and BA remains finite, although very large; (v) the 
normalization trick works nicely and gives actually the best result in a wide range of temperatures. 
The data for TAP with the normalization trick (label "TAP norm") are interrupted because at 
large f3 the iterative procedure I am using for finding the parameters {Aj} stops converging. 

The same qualitative conclusions can be reached by studying a diluted ferromagnet on a 2D 
square lattice for several different dilutions (see Figure [6]) . In particular the relative quality of the 
approximations seems to be independent on the dilution, and the BA with the normalization trick 
outperforms the other inference methods. However I notice that, while the error of BA (with and 
without the normalization trick) at the critical temperature is roughly independent on the dilution, 
the errors made by TAP and SM tend to increase when the dilution is stronger and the system 
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FIG. 7: Errors in the couplings inferred by several approximations. The model is a diluted ferromagnet on a 
3D cubic lattice of size N — 4 3 . Shown are two typical samples with different dilutions. The vertical dotted 
lines mark the loci of the ferromagnetic phase transitions in the thermodynamic limit. 

becomes more heterogenous. 

All the conclusions reached for the 2D case, perfectly apply also to the case of a diluted ferro- 
magnetic model on a 3D cubic lattice (see Figure [7]). So it is very reasonable to conclude that the 
statements made at the beginning of this Section apply to any diluted ferromagnet independently 
on the specific topology. 



B. Spin glass models 

Also for spin glass models the conclusions reached above apply very well: the relative goodness 
of the various approximations is roughly unchanged. The only major difference is that the normal- 
ization trick does not provide any more a clear improvement: its performances strongly depend on 
the disorder sample. In Figure [5] I am reporting the error on the inferred coupling for spin glass 
models on a random regular graph with fixed degree c = 4, on a 3D cubic lattice and a 2D square 
lattice. Again, willing to suggest a general purpose inference method, the choice is clearly in favor 
of the BA. Please notice that by running SuscProp it would be impossible to obtain the results 
shown in Figure El because of the limited range of convergence of such an algorithm. Indeed for 
spin glass models on a random graph SuscProp stops converging around the critical temperature 
and for spin glass models on a regular lattice it stops converging even before, well into the high 
temperature phase: for example for the a spin glass model on a 2D square lattice it converges up 
to /3bp — 0.66 25|] and on a 3D cubic lattice up to /?bp — 0.49 3l|]. In this sense, the use of the 



new formula in Eq. (|34p is really innovative. 
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FIG. 8: Errors in the couplings inferred by several approximations for a typical sample of spin glass models 
( Jij = ±/3) on different topologies: random regular graphs with fixed degree 4 (upper panels), 3D cubic and 
2D square lattices (lower panels) . The vertical dotted lines mark the loci of the spin glass phase transitions 
in the thermodynamic limit. 

C. Spin glass models with an external field 



Let me finally come to the most surprising case: frustrated models in presence of an external 
field. As shown in Figure El once more the relative level of accuracy of the 4 approximations tested 
is the same, but the is a major difference with respect to case of zero field (which has been reported 
in the upper left panel of Figure [9] for reader convenience). At high enough (3 values, both TAP 
and BA cease to have a solution to Eqs. (|31|) and (|32|) or equivalently the expressions under the 
square roots in Eqs. (|33p and (|34p become negative. To my knowledge, this fact has been never 
noticed in the past, although the TAP approximation for inferring couplings is largely used. 

In the Appendix I sketch the analytical solution to the simplest model showing this phenomenon, 
namely a system of 3 spins connected by antiferromagnetic couplings in presence of an external 
field: such an analytical solution should convince the reader that the phenomenon is not due to 
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FIG. 9: Errors in the couplings inferred by several approximations. SK model h — (upper left, just for 
comparison), SK model h = 0.1 (upper right), SK model h = 0.3 (lower left) and spin glass on a random 
regular graph with fixed degree c = 4 and h = 0.7 (lower right). The vertical dotted lines mark the loci of 
the spin glass phase transitions in the thermodynamic limit. In presence of an external field TAP and BA 
cease to have a solution for high enough j3 values. 



any numerical inaccuracy related to the complexity of the models studied here, but it can be 
mathematically proved in a very simple model. 

The absence of solutions in TAP and B A is evident in Figure [9] where the corresponding curves 
are interrupted at a (3 value that becomes smaller for larger fields. Beyond the point where BA 
stops providing inferred couplings, one should resort to other inference methods. Unfortunately at 
that point both SM and IP already give quite large errors, that keep growing fast. So in practice 
none of the methods studied in the present work is valid for inferring couplings in a frustrated 
model in presence of an external field at low enough temperatures. 
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VI. CONCLUSIONS 

The purpose of the present work is to make a detailed comparison among several approximation 
for solving the inverse Ising problem, i.e., estimating coupling and fields form magnetizations and 
correlations. 

After having explained how to derive the mean-field approximations based on the Plefka expan- 
sion (naive mean-field, TAP, Bethe approximations, etc.), I have ranked these approximations on 
the basis of how good they are in solving the direct problem (i.e., computing the correlations given 
the couplings). TAP and Bethe turned out to be in general the best approximations available. 

Secondly I have derived the new analytical formula f)34j) for inferring couplings from magne- 
tizations and correlations under the Bethe approximation. This formula allow to infer couplings 
without running the Susceptibility Propagation algorithm, thus avoiding all the serious problems 
related to the lack of convergence of such an algorithm. 

After having summarized the formulas giving the inferred couplings for the 4 approximations 
tested (independent-pair, TAP, Bethe and the small correlation expansion of Ref. [l6), I have 
introduced a trick that, normalizing the correlation matrix, improves the TAP and the Bethe 
approximations in case of models being unfrustrated or weakly frustrated. 

Finally I have presented the results of the comparison among the 4 approximations for inferring 
couplings in diluted ferromagnetic models and spin glass models. I have used several different 
topologies: fully connected graphs, regular random graphs, 3D cubic lattices and 2D square lattices. 

At the beginning of Section [V] a list of general statements about the performances of these 
approximations in solving the inverse Ising problem is given. The bottom-line suggestion is to use 
the Bethe approximation, Eq. (|34p . eventually with the normalization trick if the model is weakly 
frustrated or unfrustrated. 

In case of frustrated models with an external field (that is with non-zero magnetizations) I have 
found an important limitation for the TAP and Bethe approximations: at low enough temperatures 
these approximations stop having a solution and can not be used any more for solving the inverse 
Ising problem. This is a fundamental limitation, that take place also in very simple systems (see 
the Appendix) and that was not noticed before. 

Moreover, when the Bethe approximation stops inferring couplings, the other methods already 
have a rather large error. So, in my opinion, it is still an open problem to find an approximation 
that, using only the correlation matrix, is able to solve the inverse Ising problem in a frustrated 
model with a field at low enough temperatures. 
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Appendix A: Limits of TAP and BA inference methods for a frustrated model in a field 

In this appendix I show explicitly that the formulas derived with TAP and BA for solving the 
inverse Ising problem do not always admit a solution for the case of frustrated models. In order to 
simplify the computation I focus on the simplest model showing this problem, namely a system of 
3 spins interacting with antiferromagnetic couplings of intensity J < 0, in presence of an external 
field of intensity h, whose probability distribution is 

P(si,s 2 , s 3 ) oc exp[J(sis 2 + s 2 s 3 + s 3 si) + h(si + s 2 + S3)] • 

Thanks to the symmetries in the above measure, each spin has the same local magnetization m(J, h) 
and each pair of spins has the same correlation c( J, h). 

When using the TAP approximation for the inverse problem one has to solve the following 
equation for each coupling Jjj, 

ImimjJlj + Jij + (C -1 )jj = , 
and the above equation admit a solution only if its discriminant is non-negative: 

a tap = 1 _ 8m . m .^ C - 1 ) ij > . (Al) 

In the present case, the discriminant is the same for each coupling and it is a function of the two 
parameters J and h, that I report schematically in Figure [TOl The full curve shown in Figure [TOl 
corresponds to A TAP (J, h) = and has two asymptotes at h* = 0.966759. . . and J* = — ln(2)/4. 
It is clear that for any non-zero field h and any antiferromagnetic coupling J the inference method 
based on the TAP approximation will fail at sufficiently small temperatures (i.e., large absolute 
values of h and J). 



24 



5 



TAP 
BA 



-- 4 



- 3 



A < 



h 




- 2 



- 1 



A> 







-2 



-1.5 



-1 



-0.5 







J 



FIG. 10: Limit of validity of the TAP inference method for a system of 3 spins interacting with an antifer- 
romagnetic coupling J in an external field h. In the region where the discriminant A is negative the TAP 
inference method does not work. 

The same phenomenon happens also for the inference method based on the BA. In this case the 
discriminant that may become negative is 



In Figure [10] the dashed line corresponds to A (J, h) = and has two asymptotes at h* = 0.673689 
and along the line h = — 4 J (meaning that for this simple system of 3 spins the BA can work even 
a very low temperatures if the external field is large enough). 
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